library(readr)
library(ggplot2)
library(dplyr)
library(rms)
library(synthpop)
source('./helper functions.R')

Reading in original dataset, performing multiple imputation and writing out imputed dataset. This chunk commented out because we’ve already done this and saved the resulting file as our starting point.

Read in dataset from Brinati paper. Dataset already been imputed based on code above.

require(readr)
brinati_covid_study_v2.imputed <- read_csv("data/brinati-covid_study_v2_imputed.csv",
  col_types = cols(GENDER = col_factor(levels = c("M",
  "F")), SWAB = col_factor(levels = c("0",
  "1"))))
p<-plot_model(glm.orig.fit)
p<-plot_model(glm.orig.fit)
p<-ApplyFigureThemeLargeFontOnly(p + ylim(.1, 2.5) ) 
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

Draw ROC curve

p <- ff %>% mutate(fpr=1-Specificity) %>%
  ggplot(., aes(x=fpr, y=Sensitivity)) + geom_point(aes(frame=Cutoff)) +
    xlab('False Positive Rate (1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
Ignoring unknown aesthetics: frame
ggplotly(p)

Draw calibration

Calculate hosmer lemeshow statistics

require(performance)
hl.org <- performance_hosmer(glm.orig.fit)
hl.org
# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 5.333
           df: 8    
      p-value: 0.722
Summary: model seems to fit well.

Create smoothed cal curve from Van hoorde et al

p<-CreateSmoothedCalCurvePlot(predicted, observed)
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/smoothed-cal-curve-orig-data.png')
p

Boot strap performance in the original dataset

Generate new a synthetic dataset with a prevalence of 10%

#Assess performance of the original model on new prevalence 10 data
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(brinati.syn.prev10))), 
           observed = brinati.syn.prev10$SWAB)
          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p             U      U:Chi-sq 
 6.462721e-01  8.231360e-01 -1.737270e+00 -6.089445e-01 -1.688955e+02            NA  7.433448e-01  2.093932e+02 
          U:p             Q         Brier     Intercept         Slope          Emax           E90          Eavg 
 0.000000e+00 -1.352289e+00  2.037312e-01 -2.471404e+00  6.986839e-01  5.722161e-01  5.415682e-01  3.070988e-01 
          S:z           S:p 
 3.691032e+00  2.233457e-04 

Generate new a synthetic dataset with a prevalence of 63% – close to the original dataset

#Assess performance of the original model on new prevalence 10 data
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(brinati.syn.prev10))), 
           observed = brinati.syn.prev10$SWAB)
          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p             U      U:Chi-sq 
 7.079320e-01  8.539660e-01 -7.403297e-01 -3.535701e-01 -9.764606e+01            NA  5.720065e-01  1.615898e+02 
          U:p             Q         Brier     Intercept         Slope          Emax           E90          Eavg 
 0.000000e+00 -9.255766e-01  1.923112e-01 -2.097882e+00  7.796140e-01  4.378514e-01  4.367189e-01  2.627114e-01 
          S:z           S:p 
 3.977260e+00  6.971398e-05 

Boot strap performance of a complete re-moel of a prevalence 10% dataset

  1. Look at performance as a function of prevalence
prevalences = seq(0.01, 0.99, 0.01)
dfs_prevs <- lapply(prevalences, function (prev) brinati.syn.factory(prev = prev))
dfs.prevs.perf <- lapply(dfs_prevs, function(df) {
  assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(df))), observed = df$SWAB)
})

Calculate the different validation measures


dfs.prevs.perf.df <- data.frame(matrix(unlist(dfs.prevs.perf), nrow = length(prevalences), byrow = TRUE))
colnames(dfs.prevs.perf.df) <- names(dfs.prevs.perf[[1]])
dfs.prevs.perf.df$Prevalence = prevalences
dfs.prevs.perf.df <- gather(dfs.prevs.perf.df, Measure, Value, Dxy:`S:p`, factor_key = TRUE)

Plot the distribution of validation measures

measures = c('C (ROC)', 'Brier', 'Intercept', 'Slope')
ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Measure, y=Value,  color=Measure)) + geom_boxplot() +
  geom_point(data=data.frame(Measure = measures, Value = c(1, 0, 0, 1)), size=3, color='blue') + theme_bw()


ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Prevalence, y=Value,  color=Measure, group=Measure)) + geom_line() + theme_bw()

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7ciBtZXNzYWdlPUZBTFNFfQpsaWJyYXJ5KHJlYWRyKQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkocm1zKQpsaWJyYXJ5KHN5bnRocG9wKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KHBsb3RseSkKbGlicmFyeShzalBsb3QpCnNvdXJjZSgnLi9oZWxwZXIgZnVuY3Rpb25zLlInKQpgYGAKClJlYWRpbmcgaW4gb3JpZ2luYWwgZGF0YXNldCwgcGVyZm9ybWluZyBtdWx0aXBsZSBpbXB1dGF0aW9uIGFuZCB3cml0aW5nIG91dCBpbXB1dGVkIGRhdGFzZXQuIFRoaXMgY2h1bmsgY29tbWVudGVkIG91dCBiZWNhdXNlIHdlJ3ZlIGFscmVhZHkgZG9uZSB0aGlzIGFuZCBzYXZlZCB0aGUgcmVzdWx0aW5nIGZpbGUgYXMgb3VyIHN0YXJ0aW5nIHBvaW50LgpgYGB7ciBlY2hvPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQojIAojIGJyaW5hdGlfY292aWRfc3R1ZHlfdjIgPC0gcmVhZF9jc3YoImRhdGEvYnJpbmF0aS1jb3ZpZF9zdHVkeV92Mi5jc3YiLAojIGNvbF90eXBlcyA9IGNvbHMoR0VOREVSID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCJNIiwKIyAiRiIpKSwgU1dBQiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMCIsCiMgIjEiKSkpKQojIGJyaW5hdGlfY292aWRfc3R1ZHlfdjIgPC0gYnJpbmF0aV9jb3ZpZF9zdHVkeV92MiAlPiUgbXV0YXRlX2lmKGlzLm51bWVyaWMsIC5mdW5zID0gZnVuY3Rpb24oeCkge3grMC4wMDF9KQojIGxpYnJhcnkobWkpCiMgYnJhaW5hdGkubWkgPC0gbWlzc2luZ19kYXRhLmZyYW1lKGFzLmRhdGEuZnJhbWUoYnJpbmF0aV9jb3ZpZF9zdHVkeV92MiksIGZhdm9yX3Bvc2l0aXZlID0gVFJVRSkKIyAKIyAKIyBicmFpbmF0aS5taTwtY2hhbmdlKGJyYWluYXRpLm1pLCB5PWMoJ0x5bXBob2N5dGVzJywgJ0Jhc29waGlscycsICdNb25vY3l0ZXMnLCdFb3Npbm9waGlscycsICdCYXNvcGhpbHMnKSwgd2hhdD0ndHlwZScsIHRvPSdwb3NpdGl2ZS1jb250aW51b3VzJykKIyBzaG93KGJyYWluYXRpLm1pKQojIGltYWdlKGJyYWluYXRpLm1pKQojIAojIG9wdGlvbnMobWMuY29yZXMgPSA0KQojIGltcHV0YXRpb25zIDwtIG1pKGJyYWluYXRpLm1pLCBuLml0ZXIgPSA5MCwgbi5jaGFpbnMgPSA0LCBtYXgubWludXRlcyA9IDIwKSAKIyBzaG93KGltcHV0YXRpb25zKQojIHJvdW5kKG1pcHBseShpbXB1dGF0aW9ucywgbWVhbiwgdG8ubWF0cml4ID0gVFJVRSksIDMpCiMgUmhhdHMoaW1wdXRhdGlvbnMpCiMgZGZzIDwtIGNvbXBsZXRlKGltcHV0YXRpb25zLCBtPTIpCiMgYnJpbmF0aV9jb3ZpZF92Mi5pbXB1dGVkIDwtIGRmc1tbMV1dWywxOjE2XQojIHdyaXRlLmNzdihicmluYXRpX2NvdmlkX3YyLmltcHV0ZWQsIGZpbGU9Jy4vZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2Jywgcm93Lm5hbWVzID0gRkFMU0UpCmBgYAoKUmVhZCBpbiBkYXRhc2V0IGZyb20gQnJpbmF0aSBwYXBlci4gRGF0YXNldCBhbHJlYWR5IGJlZW4gaW1wdXRlZCBiYXNlZCBvbiBjb2RlIGFib3ZlLiAKYGBge3J9CnJlcXVpcmUocmVhZHIpCmJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCA8LSByZWFkX2NzdigiZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2IiwKICBjb2xfdHlwZXMgPSBjb2xzKEdFTkRFUiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiTSIsCiAgIkYiKSksIFNXQUIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjAiLAogICIxIikpKSkKYGBgCgoKYGBge3J9CgpnbG0ub3JpZy5maXQ8LWdsbShTV0FCIH4gLiwgYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkLCBmYW1pbHk9J2Jpbm9taWFsJykKZ2xtLm9yaWcuZml0CndyaXRlLmNzdihzdW1tYXJ5KGdsbS5vcmlnLmZpdCkkY29lZiwgZmlsZT0nb3V0cHV0L09SLW9yaWdpbmFsLW1vZGVsLmNzdicsIHJvdy5uYW1lcyA9IEZBTFNFKQpwcmVkaWN0ZWQgPSBwbG9naXMocHJlZGljdChnbG0ub3JpZy5maXQpKQpvYnNlcnZlZCA9IGJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCRTV0FCCmFzc2Vzc1BlcmYocHJlZGljdGVkLCBvYnNlcnZlZCkKcDwtcGxvdF9tb2RlbChnbG0ub3JpZy5maXQpCnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHAgKyB5bGltKC4xLCAyLjUpICkgCiAgCnAKZ2dzYXZlKGZpbGVuYW1lID0gJ2ZpZ3MvT1Itb3JpZ2luYWwtbW9kZWwucG5nJywgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkKYGBgCgojIERyYXcgUk9DIGN1cnZlCmBgYHtyfQpjbXMuY3V0b2ZmcyA8LSBsYXBwbHkoc2VxKDAuMDEsMC45OSwwLjAxKSwgZnVuY3Rpb24oY3V0b2ZmKSB7CiAgcmV0PWNvbmZ1c2lvbk1hdHJpeChwcmVkaWN0ZWQsIG9ic2VydmVkLCBjdXRvZmY9Y3V0b2ZmKQogIHJldCRjdXRvZmYgPSBjdXRvZmYKICByZXQKICB9KQpmZjwtZGF0YS5mcmFtZSh0KHNhcHBseShjbXMuY3V0b2ZmcywgZnVuY3Rpb24oY2YpIGMoIkN1dG9mZiI9Y2YkY3V0b2ZmLCAiU2Vuc2l0aXZ0eSI9Y2Ykc2VucywgIlNwZWNpZmljaXR5Ij1jZiRzcGVjKSkpKQpjb2xuYW1lcyhmZikgPC0gYygnQ3V0b2ZmJywgJ1NlbnNpdGl2aXR5JywnU3BlY2lmaWNpdHknKQoKCmN1dG9mZnNfdG9fcGxvdCA8LSBjKDAuMDksIDAuMywgMC42LCAwLjkpCnAgPC0gQXBwbHlGaWd1cmVUaGVtZShmZiAlPiUgbXV0YXRlKGZwcj0xLVNwZWNpZmljaXR5KSAlPiUgZmlsdGVyKEN1dG9mZiAlaW4lIGN1dG9mZnNfdG9fcGxvdCkgJT4lCiAgZ2dwbG90KC4sIGFlcyh4PWZwciwgeT1TZW5zaXRpdml0eSkpICsgCiAgICBnZW9tX3BvaW50KHNpemU9MywgY29sb3I9J2RhcmtibHVlJykgKyB4bGltKDAsMSkgKyB5bGltKDAsMSkrCiAgICB4bGFiKCdGYWxzZSBQb3NpdGl2ZSBSYXRlXG4oMS1TcGVjaWZpY2l0eSknKSArIHlsYWIoJ1RydWUgUG9zaXRpdmUgUmF0ZVxuKFNlbnNpdGl2aXR5KScpCiAgKQpnZ3NhdmUoZmlsZW5hbWUgPSAnZmlncy9leGFtcGxlLXNlbnMtc3BlYy1vbi1yb2MucG5nJywgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkKCmZvciAoY3V0b2ZmIGluIGxhcHBseSgxOmxlbmd0aChjdXRvZmZzX3RvX3Bsb3QpLCBmdW5jdGlvbihpKSBjdXRvZmZzX3RvX3Bsb3RbMTppXSkpIHsKICBwIDwtIEFwcGx5RmlndXJlVGhlbWUoZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lIGZpbHRlcihDdXRvZmYgJWluJSBjdXRvZmYpICU+JQogICAgZ2dwbG90KC4sIGFlcyh4PWZwciwgeT1TZW5zaXRpdml0eSkpICsgCiAgICAgIGdlb21fcG9pbnQoc2l6ZT0zLCBjb2xvcj0nZGFya2JsdWUnKSArIHhsaW0oMCwxKSArIHlsaW0oMCwxKSsKICAgICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQogICkKICBnZ3NhdmUoZmlsZW5hbWUgPSBwYXN0ZTAoJ2ZpZ3MvZXhhbXBsZS1zZW5zLXNwZWMtb24tcm9jLScsIGN1dG9mZltsZW5ndGgoY3V0b2ZmKV0sICAnLnBuZycpLCAKICAgICAgICAgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkKfQpwIDwtIEFwcGx5RmlndXJlVGhlbWUoCiAgZ2dwbG90KGZmLCBhZXMoeD0xLVNwZWNpZmljaXR5LCB5PVNlbnNpdGl2aXR5KSkgKyBnZW9tX2xpbmUoKSsgCiAgICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IDAsIHltYXg9U2Vuc2l0aXZpdHkpLCBjb2xvcj1OQSwgZmlsbD0nYmx1ZScsYWxwaGE9MC4zKSArIAogICAgeGxpbSgwLDEpICsgeWxpbSgwLDEpKwogICAgICB4bGFiKCdGYWxzZSBQb3NpdGl2ZSBSYXRlXG4oMS1TcGVjaWZpY2l0eSknKSArIHlsYWIoJ1RydWUgUG9zaXRpdmUgUmF0ZVxuKFNlbnNpdGl2aXR5KScpIAopCmdnc2F2ZShmaWxlbmFtZSA9ICdmaWdzL2V4YW1wbGUtcm9jLnBuZycsIAogICAgICAgICBwbG90PXAsIHdpZHRoPTYsIGhlaWdodD02LCB1bml0cz0naW4nLCBkcGk9MzAwKQpwIDwtIGZmICU+JSBtdXRhdGUoZnByPTEtU3BlY2lmaWNpdHkpICU+JQogIGdncGxvdCguLCBhZXMoeD1mcHIsIHk9U2Vuc2l0aXZpdHkpKSArIGdlb21fcG9pbnQoYWVzKGZyYW1lPUN1dG9mZikpICsKICAgIHhsYWIoJ0ZhbHNlIFBvc2l0aXZlIFJhdGUgKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQogIAogIApnZ3Bsb3RseShwKQpgYGAKCiMgRHJhdyBjYWxpYnJhdGlvbgpgYGB7cn0KY2FsX2JyZWFrcyA9IHNlcSgwLCAxLCAuMSkKY2FsX2RlY2lsZXMgPC0gZGF0YS5mcmFtZShQcmVkaWN0ZWQ9cHJlZGljdGVkLCBPYnNlcnZlZCA9IGFzLmludGVnZXIoYXMuY2hhcmFjdGVyKG9ic2VydmVkKSkpICU+JQogIG11dGF0ZShCaW5zID0gY3V0KFByZWRpY3RlZCwgYnJlYWtzPWNhbF9icmVha3MsIGluY2x1ZGUubG93ZXN0ID0gVFJVRSkpCmNhbF9kZWNpbGVzX3N1bW1hcnkgPC0gY2FsX2RlY2lsZXMgJT4lIGdyb3VwX2J5KEJpbnMpICU+JSBzdW1tYXJpc2Uobj1uKCksIEF2ZXJhZ2VfT2JzZXJ2ZWQgPSBtZWFuKE9ic2VydmVkKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBBdmVyYWdlX1ByZWRpY3RlZD1tZWFuKFByZWRpY3RlZCkpIAoKCmZvciAoY3V0b2ZmIGluIGxhcHBseSgxOm5yb3coY2FsX2RlY2lsZXNfc3VtbWFyeSksIGZ1bmN0aW9uKGkpIGNhbF9kZWNpbGVzX3N1bW1hcnkkQmluc1sxOmldKSkgewogIHA8LUFwcGx5RmlndXJlVGhlbWVDYWxDdXJ2ZVBvaW50cygKICAgIGdncGxvdChjYWxfZGVjaWxlc19zdW1tYXJ5ICU+JSBmaWx0ZXIoQmlucyAlaW4lIGN1dG9mZiksIGFlcyh4PUF2ZXJhZ2VfUHJlZGljdGVkLCB5PUF2ZXJhZ2VfT2JzZXJ2ZWQpKSApCiAgcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkKICBTYXZlU3RkU3F1YXJlRmlndXJlKHAsIGZpbGVuYW1lID0gcGFzdGUwKCdmaWdzL2V4YW1wbGUtY2FsLWN1cnZlLXBvaW50cy0nLCBjdXRvZmZbbGVuZ3RoKGN1dG9mZildLCAgJy5wbmcnKSkKfQoKcDwtQXBwbHlGaWd1cmVUaGVtZUNhbEN1cnZlUG9pbnRzKGdncGxvdChjYWxfZGVjaWxlc19zdW1tYXJ5ICwgYWVzKHg9QXZlcmFnZV9QcmVkaWN0ZWQsIHk9QXZlcmFnZV9PYnNlcnZlZCkpICkKcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkKClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3MvY2FsLWN1cnZlLWRlY2lsZXMtYWxsLnBuZycpCgpwPC1wICsgZ2VvbV9zbW9vdGgobWV0aG9kPSdsbScsIHNlPUZBTFNFKQpTYXZlU3RkU3F1YXJlRmlndXJlKHAsICdmaWdzL2NhbC1jdXJ2ZS1kZWNpbGVzLWFsbC13aXRoLWJlc3RmaXQucG5nJykKcAogICAgICAgICAgICAKYGBgCgpDYWxjdWxhdGUgaG9zbWVyIGxlbWVzaG93IHN0YXRpc3RpY3MKYGBge3J9CnJlcXVpcmUocGVyZm9ybWFuY2UpCmhsLm9yZyA8LSBwZXJmb3JtYW5jZV9ob3NtZXIoZ2xtLm9yaWcuZml0KQpobC5vcmcKYGBgCgpDcmVhdGUgc21vb3RoZWQgY2FsIGN1cnZlIGZyb20gVmFuIGhvb3JkZSBldCBhbApgYGB7cn0KcDwtQ3JlYXRlU21vb3RoZWRDYWxDdXJ2ZVBsb3QocHJlZGljdGVkLCBvYnNlcnZlZCkKcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkKU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9zbW9vdGhlZC1jYWwtY3VydmUtb3JpZy1kYXRhLnBuZycpCnAKYGBgCgpCb290IHN0cmFwIHBlcmZvcm1hbmNlIGluIHRoZSBvcmlnaW5hbCBkYXRhc2V0CmBgYHtyIHdhcm5pbmc9RkFMU0UsIGVjaG89RkFMU0V9CmRmIDwtIGFzLmRhdGEuZnJhbWUoYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkKQpuQm9vdCA9IDEwMApzZXQuc2VlZCgxMzQyMzQ5KQp0cmFpbklkeEJvb3QgPC0gbGFwcGx5KDE6bkJvb3QsIGZ1bmN0aW9uKGkpIHNhbXBsZSgxOm5yb3coZGYpLCBzaXplPW5yb3coZGYpLCByZXBsYWNlPVRSVUUpKQoKYm9vdFBlcmYub3JpZ2RhdGEgPC0gbGFwcGx5KHRyYWluSWR4Qm9vdCwgCiAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih0cmFpbklkeCkgCiAgICAgICAgICAgICAgICAgICAgIGFzc2Vzc1NpbmdsZVRyYWluVGVzdCh0cmFpbklkeCwgKDE6bnJvdyhkZikpWy10cmFpbklkeF0sIGRmLCBnbG0ubW9kZWwgPSBTV0FCfiAuLCBvdXRjb21lLnZhci5uYW1lID0gJ1NXQUInKQogICAgICAgICAgICAgICAgICAgKSAKCm1lYXN1cmVzID0gYygnQyAoUk9DKScsICdCcmllcicsICdJbnRlcmNlcHQnLCAnU2xvcGUnKQpicmluYXRpLm9yaWcuYm9vdC5wZXJmIDwtIGRhdGEuZnJhbWUobWF0cml4KAogIHVubGlzdChsYXBwbHkoYm9vdFBlcmYub3JpZ2RhdGEsIGZ1bmN0aW9uKGJvb3QpIGJvb3QkdGVzdC5wZXJmKSksIAogIG5yb3cgPSBuQm9vdCwgYnlyb3cgPSBUUlVFKSkKY29sbmFtZXMoYnJpbmF0aS5vcmlnLmJvb3QucGVyZikgPC0gbmFtZXMoYm9vdFBlcmYub3JpZ2RhdGFbWzFdXSR0ZXN0LnBlcmYpCgpicmluYXRpLm9yaWcuYm9vdC5wZXJmIDwtIGdhdGhlcihicmluYXRpLm9yaWcuYm9vdC5wZXJmLCBNZWFzdXJlLCBWYWx1ZSwgRHh5OmBTOnBgLCBmYWN0b3Jfa2V5ID0gVFJVRSkKcDwtZ2dwbG90KGJyaW5hdGkub3JpZy5ib290LnBlcmYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgKwogIGdlb21fcG9pbnQoZGF0YT1kYXRhLmZyYW1lKE1lYXN1cmUgPSBtZWFzdXJlcywgVmFsdWUgPSBjKDEsIDAsIDAsIDEpKSwgc2l6ZT0zLCBjb2xvcj0nYmx1ZScpICsgdGhlbWVfYncoKQpwPC1BcHBseUZpZ3VyZVRoZW1lKHApClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3Mvb3JpZy1tb2RlbC1ib290c3RyYXAtZXZhbHVhdGlvbi5wbmcnKQpwCmBgYAoKCkdlbmVyYXRlIG5ldyBhIHN5bnRoZXRpYyBkYXRhc2V0IHdpdGggYSBwcmV2YWxlbmNlIG9mIDEwJQpgYGB7cn0KCmJyaW5hdGkuc3luLmZhY3RvcnkgPC0gZ2VuZXJhdGVTeW50aGV0aWNEYXRhRmFjdG9yeShicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQsIG1ldGhvZCA9ICdzeW4nKQpicmluYXRpLnN5bi5wcmV2MTAgPC0gYnJpbmF0aS5zeW4uZmFjdG9yeShwcmV2ID0gMC4xMCwgc2FtcFNpemUgPSAyNzkpCgoKI0Fzc2VzcyBwZXJmb3JtYW5jZSBvZiB0aGUgb3JpZ2luYWwgbW9kZWwgb24gbmV3IHByZXZhbGVuY2UgMTAgZGF0YQphc3Nlc3NQZXJmKHByZWRpY3RlZCA9IHBsb2dpcyhwcmVkaWN0KGdsbS5vcmlnLmZpdCwgbmV3ZGF0YT1hcy5kYXRhLmZyYW1lKGJyaW5hdGkuc3luLnByZXYxMCkpKSwgCiAgICAgICAgICAgb2JzZXJ2ZWQgPSBicmluYXRpLnN5bi5wcmV2MTAkU1dBQikKCiMgcmVmaXQgdGhlIG1vZGVsCmdsbS5zeW4uZml0PC1nbG0oU1dBQiB+IC4sIGJyaW5hdGkuc3luLnByZXYxMCwgZmFtaWx5PSdiaW5vbWlhbCcpCmdsbS5zeW4uZml0CmBgYApHZW5lcmF0ZSBuZXcgYSBzeW50aGV0aWMgZGF0YXNldCB3aXRoIGEgcHJldmFsZW5jZSBvZiA2MyUgLS0gY2xvc2UgdG8gdGhlIG9yaWdpbmFsIGRhdGFzZXQKYGBge3J9CgpicmluYXRpLnN5bi5mYWN0b3J5IDwtIGdlbmVyYXRlU3ludGhldGljRGF0YUZhY3RvcnkoYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkLCBtZXRob2QgPSAnYm9vdCcpCmJyaW5hdGkuc3luLnByZXYxMCA8LSBicmluYXRpLnN5bi5mYWN0b3J5KHByZXYgPSAwLjEsIHNhbXBTaXplID0gMjc5KSAlPiUgc2VsZWN0KC13ZWlnaHQpCgoKI0Fzc2VzcyBwZXJmb3JtYW5jZSBvZiB0aGUgb3JpZ2luYWwgbW9kZWwgb24gbmV3IHByZXZhbGVuY2UgMTAgZGF0YQphc3Nlc3NQZXJmKHByZWRpY3RlZCA9IHBsb2dpcyhwcmVkaWN0KGdsbS5vcmlnLmZpdCwgbmV3ZGF0YT1hcy5kYXRhLmZyYW1lKGJyaW5hdGkuc3luLnByZXYxMCkpKSwgCiAgICAgICAgICAgb2JzZXJ2ZWQgPSBicmluYXRpLnN5bi5wcmV2MTAkU1dBQikKCiMgcmVmaXQgdGhlIG1vZGVsCmdsbS5zeW4uZml0PC1nbG0oU1dBQiB+IC4sIGJyaW5hdGkuc3luLnByZXYxMCwgZmFtaWx5PSdiaW5vbWlhbCcpCmdsbS5zeW4uZml0CmBgYAoKQm9vdCBzdHJhcCBwZXJmb3JtYW5jZSBvZiBhIGNvbXBsZXRlIHJlLW1vZWwgb2YgYSBwcmV2YWxlbmNlIDEwJSBkYXRhc2V0CmBgYHtyIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmdsbS5tb2RlbCA9IFNXQUIgfiAuIApuQm9vdCA9IDEwMApkZiA8LSBhcy5kYXRhLmZyYW1lKGJyaW5hdGkuc3luLnByZXYxMCkKc2V0LnNlZWQoMTM0MjM0OSkKdHJhaW5JZHhCb290IDwtIGxhcHBseSgxOm5Cb290LCBmdW5jdGlvbihpKSBzYW1wbGUoMTpucm93KGRmKSwgc2l6ZT1ucm93KGRmKSwgcmVwbGFjZT1UUlVFKSkKCmJvb3RQZXJmLnByZXYxMCA8LSBsYXBwbHkodHJhaW5JZHhCb290LCAKICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHRyYWluSWR4KSAKICAgICAgICAgICAgICAgICAgICAgYXNzZXNzU2luZ2xlVHJhaW5UZXN0KHRyYWluSWR4LCAoMTpucm93KGRmKSlbLXRyYWluSWR4XSwgZGYsIGdsbS5tb2RlbCwgb3V0Y29tZS52YXIubmFtZSA9ICdTV0FCJykKICAgICAgICAgICAgICAgICAgICkgCmF1Yy5ib290IDwtIHNhcHBseShib290UGVyZi5wcmV2MTAsIGZ1bmN0aW9uKGJvb3QpIGJvb3QkdGVzdC5wZXJmWydDIChST0MpJ10pCnNwcmludGYoIkFVQzogJTAuMmYgKCUwLjJmLCAlMC4yZikiLCBtZWRpYW4oYXVjLmJvb3QpLCBxdWFudGlsZShhdWMuYm9vdCwgcHJvYnM9YygwLjAyNSkpLCBxdWFudGlsZShhdWMuYm9vdCwgcHJvYnM9YygwLjk3NSkpKQpgYGAKCklJLiBMb29rIGF0IHBlcmZvcm1hbmNlIGFzIGEgZnVuY3Rpb24gb2YgcHJldmFsZW5jZQpgYGB7cn0KcHJldmFsZW5jZXMgPSBzZXEoMC4wMSwgMC45OSwgMC4wMSkKZGZzX3ByZXZzIDwtIGxhcHBseShwcmV2YWxlbmNlcywgZnVuY3Rpb24gKHByZXYpIGJyaW5hdGkuc3luLmZhY3RvcnkocHJldiA9IHByZXYpKQpkZnMucHJldnMucGVyZiA8LSBsYXBwbHkoZGZzX3ByZXZzLCBmdW5jdGlvbihkZikgewogIGFzc2Vzc1BlcmYocHJlZGljdGVkID0gcGxvZ2lzKHByZWRpY3QoZ2xtLm9yaWcuZml0LCBuZXdkYXRhPWFzLmRhdGEuZnJhbWUoZGYpKSksIG9ic2VydmVkID0gZGYkU1dBQikKfSkKYGBgCgpDYWxjdWxhdGUgdGhlIGRpZmZlcmVudCB2YWxpZGF0aW9uIG1lYXN1cmVzCmBgYHtyfQoKZGZzLnByZXZzLnBlcmYuZGYgPC0gZGF0YS5mcmFtZShtYXRyaXgodW5saXN0KGRmcy5wcmV2cy5wZXJmKSwgbnJvdyA9IGxlbmd0aChwcmV2YWxlbmNlcyksIGJ5cm93ID0gVFJVRSkpCmNvbG5hbWVzKGRmcy5wcmV2cy5wZXJmLmRmKSA8LSBuYW1lcyhkZnMucHJldnMucGVyZltbMV1dKQpkZnMucHJldnMucGVyZi5kZiRQcmV2YWxlbmNlID0gcHJldmFsZW5jZXMKZGZzLnByZXZzLnBlcmYuZGYgPC0gZ2F0aGVyKGRmcy5wcmV2cy5wZXJmLmRmLCBNZWFzdXJlLCBWYWx1ZSwgRHh5OmBTOnBgLCBmYWN0b3Jfa2V5ID0gVFJVRSkKYGBgCgpQbG90IHRoZSBkaXN0cmlidXRpb24gb2YgdmFsaWRhdGlvbiBtZWFzdXJlcwpgYGB7cn0KbWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpCmdncGxvdChkZnMucHJldnMucGVyZi5kZiAlPiUgZmlsdGVyKE1lYXN1cmUgJWluJSBtZWFzdXJlcyksIGFlcyh4PU1lYXN1cmUsIHk9VmFsdWUsICBjb2xvcj1NZWFzdXJlKSkgKyBnZW9tX2JveHBsb3QoKSArCiAgZ2VvbV9wb2ludChkYXRhPWRhdGEuZnJhbWUoTWVhc3VyZSA9IG1lYXN1cmVzLCBWYWx1ZSA9IGMoMSwgMCwgMCwgMSkpLCBzaXplPTMsIGNvbG9yPSdibHVlJykgKyB0aGVtZV9idygpCgpnZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1QcmV2YWxlbmNlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSwgZ3JvdXA9TWVhc3VyZSkpICsgZ2VvbV9saW5lKCkgKyB0aGVtZV9idygpCmBgYAoK